In this exercise we’ll practice joining climate data tables from Cal-Adapt to other tables so we can create indices and map results. Our exercise will focus on census tracts in Kern County, and the end goal will be to create maps with a composite index of census tract cumulative impact score from CalEnviroScreen (a measure of exposure to pollution) and the anticipatd increase in temperature from historic condition to the end-of-century under RCP8.5.
First we load the packages we need for this exercise, which include caladaptr, dplyr, and sf:
library(caladaptr)
library(sf)
library(dplyr)
library(lubridate)
library(units)
library(ggplot2)
library(tidyr)
library(scales)
library(tmap)
Start by downloading all the census tracts:
tracts_sf <- ca_aoipreset_geom("censustracts")
Reading layer `censustracts' from data source `C:\Users\Andy\AppData\Local\R\cache\R\caladaptr\censustracts.gpkg' using driver `GPKG'
Simple feature collection with 8034 features and 63 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -13850030 ymin: 3833633 xmax: -12705030 ymax: 5162403
projected CRS: WGS 84 / Pseudo-Mercator
head(tracts_sf)
Simple feature collection with 6 features and 63 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -13610980 ymin: 4556079 xmax: -13604620 ymax: 4563255
projected CRS: WGS 84 / Pseudo-Mercator
tract pop2010 ciscore ciscorep ces_3_0_percentile_range ozone ozonep pm pmp
1 6001400100 2937 3 2 1-5% (lowest scores) 0 8 8 31
2 6001400200 1974 1 0 1-5% (lowest scores) 0 8 8 31
3 6001400300 4865 12 18 15-20% 0 8 8 31
4 6001400400 3703 10 15 15-20% 0 8 8 31
5 6001400500 3517 18 35 30-35% 0 8 8 31
6 6001400600 1571 18 34 30-35% 0 8 8 31
diesel dieselp drink drinkp pest pestp rseihaz rseihazp traffic trafficp cleanups
1 27 81 70 4 0 0 485 50 929 64 4
2 42 94 70 4 0 0 442 49 1392 81 0
3 42 94 70 4 0 0 426 48 1207 75 0
4 42 94 70 4 0 0 444 49 1153 74 0
5 42 94 70 4 0 0 448 49 668 48 3
6 42 94 70 4 0 0 426 48 1743 88 1
cleanupsp gwthreats gwthreatsp haz hazp iwb iwbp swis swisp pollution pollutionscore
1 42 4 27 6 98 2 29 1 33 37 4
2 0 20 72 0 0 0 0 0 0 30 3
3 13 35 85 0 0 0 0 0 0 31 3
4 0 38 87 0 0 0 0 0 0 30 3
5 35 21 73 0 0 0 0 0 0 29 3
6 24 20 72 0 0 0 0 0 0 32 3
pollutionp asthma asthmap lbw lbwp cvd cvdp edu edup ling lingp pov povp unemp unempp
1 40 18 6 2 9 3 1 1 4 3 29 7 4 NA NA
2 20 24 15 1 0 3 5 1 3 0 0 11 9 2.2 1
3 23 37 38 4 29 4 8 5 20 8 55 14 15 8.8 46
4 21 60 70 3 26 4 12 6 24 2 17 16 20 3.3 4
5 18 110 95 5 68 7 44 3 11 2 17 25 36 6.5 25
6 25 122 97 1 1 8 55 7 31 1 8 25 37 14.9 84
housingb housingbp popchar popcharscore popcharp children_under10_pct
1 8 8 8 0 1 9
2 4 1 4 0 0 12
3 16 40 30 3 20 11
4 15 36 27 2 17 11
5 19 59 49 5 49 9
6 13 26 44 4 41 7
pop_11_64_years_pct elderly_over_65_pct hispanic_pct white_pct african_american_pct
1 69 21 4 70 4
2 71 16 7 78 1
3 78 10 8 66 10
4 79 9 9 65 12
5 82 8 9 50 26
6 82 9 8 42 39
native_american_pct asian_american_pct other_pct cidecile civigintile id
1 0 15 4 1 1 1
2 0 7 5 1 1 2
3 0 8 5 2 4 3
4 0 7 5 2 4 4
5 0 6 6 4 7 5
6 0 5 4 4 7 6
geom
1 MULTIPOLYGON (((-13608465 4...
2 MULTIPOLYGON (((-13609690 4...
3 MULTIPOLYGON (((-13610548 4...
4 MULTIPOLYGON (((-13610124 4...
5 MULTIPOLYGON (((-13610979 4...
6 MULTIPOLYGON (((-13610819 4...
We can pull out just the tracts for Kern County by decomposing the tract id. All census tracts in Kern County will have an idea number that starts with “6029” (06 is for California, 029 is Kern County’s FIPs code).
kern_tracts_sf <- tracts_sf %>%
filter(tract %>% as.character() %>% stringr::str_starts("6029"))
View the attribute table. You’ll notice the census tracts have values from CalEnviroScreen, including the cumulative impact score (ciscore), which we’ll come back to you later.
head(kern_tracts_sf)
Simple feature collection with 6 features and 63 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -13260680 ymin: 4217755 xmax: -13248100 ymax: 4228472
projected CRS: WGS 84 / Pseudo-Mercator
tract pop2010 ciscore ciscorep ces_3_0_percentile_range ozone ozonep pm pmp
1 6029000101 12240 47 86 85-90% 0 91 18 100
2 6029000102 3055 33 66 65-70% 0 91 18 100
3 6029000200 7644 50 89 85-90% 0 91 18 100
4 6029000300 4217 47 86 85-90% 0 91 18 100
5 6029000400 4319 57 95 95-100% (highest scores) 0 91 18 100
6 6029000503 6347 15 28 25-30% 0 91 18 100
diesel dieselp drink drinkp pest pestp rseihaz rseihazp traffic trafficp cleanups
1 21 66 485 52 137 78 87 24 512 32 31
2 21 65 288 31 0 0 21 12 209 6 11
3 21 67 288 31 0 0 53 20 528 34 21
4 21 66 288 31 0 0 49 19 406 22 17
5 49 96 786 89 0 0 89 24 761 54 14
6 26 79 691 82 1352 91 496 51 984 67 1
cleanupsp gwthreats gwthreatsp haz hazp iwb iwbp swis swisp pollution pollutionscore
1 91 80 96 11 100 0 0 2 50 64 7
2 65 29 81 0 43 0 0 0 0 41 5
3 85 24 77 3 95 0 0 0 0 49 6
4 81 8 43 2 93 0 0 0 0 45 5
5 75 4 29 2 93 0 0 0 0 56 6
6 18 0 4 0 26 0 0 0 20 57 7
pollutionp asthma asthmap lbw lbwp cvd cvdp edu edup ling lingp pov povp unemp unempp
1 96 80 86 5 61 11 87 11 44 1 10 30 46 8.7 45
2 50 80 86 5 64 11 87 15 53 0 0 49 72 17.4 91
3 72 80 86 6 86 11 87 22 67 2 14 75 97 16.7 90
4 60 80 86 7 96 11 87 26 71 3 23 68 92 23.6 98
5 86 79 85 6 80 11 84 24 70 4 35 76 97 21.3 97
6 88 33 31 2 4 7 46 10 42 1 12 11 9 5.3 15
housingb housingbp popchar popcharscore popcharp children_under10_pct
1 16 44 57 6 62 17
2 12 24 63 6 70 15
3 33 93 79 8 92 16
4 24 74 80 8 93 16
5 32 91 80 8 93 17
6 6 4 21 2 10 14
pop_11_64_years_pct elderly_over_65_pct hispanic_pct white_pct african_american_pct
1 74 7 19 74 0
2 69 15 11 82 0
3 75 7 21 73 0
4 71 12 21 73 0
5 73 9 21 72 1
6 76 8 17 74 0
native_american_pct asian_american_pct other_pct cidecile civigintile id
1 1 2 2 9 18 965
2 1 0 3 7 14 966
3 1 0 3 9 18 967
4 1 0 2 9 18 968
5 1 0 3 10 20 969
6 0 3 2 3 6 970
geom
1 MULTIPOLYGON (((-13256411 4...
2 MULTIPOLYGON (((-13249411 4...
3 MULTIPOLYGON (((-13253779 4...
4 MULTIPOLYGON (((-13252464 4...
5 MULTIPOLYGON (((-13253750 4...
6 MULTIPOLYGON (((-13260684 4...
For mapping, we’ll also import the county boundary for Kern:
kern_bnd_sf <- ca_aoipreset_geom("counties") %>%
filter(fips == "06029")
Reading layer `counties' from data source `C:\Users\Andy\AppData\Local\R\cache\R\caladaptr\counties.gpkg' using driver `GPKG'
Simple feature collection with 87 features and 54 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -13871160 ymin: 3833648 xmax: -12625080 ymax: 5416187
projected CRS: WGS 84 / Pseudo-Mercator
Plot them to make sure everything looks right:
library(tmap)
tmap_mode("plot")
tm_shape(kern_tracts_sf) +
tm_polygons() +
tm_shape(kern_bnd_sf) +
tm_borders(col = "red", lty = 2) +
tm_layout(main.title = "Census Tracts in Kern County, CA")
Now we can create two API requests for the historic values and end of century values:
## Create an API request for historic values of tasmax, using a 32 model ensemble
kern_trcts_hist_cap <- ca_loc_aoipreset(type = "censustracts",
idfld = "tract",
idval = kern_tracts_sf %>% pull(tract)) %>%
ca_cvar("tasmax") %>%
ca_period("year") %>%
ca_gcm("ens32max") %>%
ca_scenario("historical") %>%
ca_years(start = 1985, end = 2005) %>%
ca_options(spatial_ag = "mean")
kern_trcts_hist_cap
Cal-Adapt API Request
Location(s):
AOI Preset: censustracts
tract(s): 6029000101, 6029000102, 6029000200, ...
Variable(s): tasmax
Temporal aggregration period(s): year
GCM(s): ens32max
Scenario(s): historical
Slug(s): NA
Dates: 1985-01-01 to 2005-12-31
Options:
spatial ag: mean
temporal ag (add'l): NA
plot(kern_trcts_hist_cap)
Fetch the historic data:
kern_trcts_hist_tbl <- kern_trcts_hist_cap %>% ca_getvals_tbl(quiet = TRUE)
## backup: kern_trcts_hist_tbl <- readRDS("data/kern_trcts_hist_tbl.rds")
dim(kern_trcts_hist_tbl)
[1] 3171 8
head(kern_trcts_hist_tbl)
Next compute the average temperature for each tract.
kern_trcts_mean_temp_hist <- kern_trcts_hist_tbl %>%
group_by(tract) %>%
summarise(mean_temp_hist = mean(val))
`summarise()` ungrouping output (override with `.groups` argument)
kern_trcts_mean_temp_hist
Do the same for the end-of-century period with RCP85.
kern_trcts_prj_cap <- ca_loc_aoipreset(type = "censustracts",
idfld = "tract",
idval = kern_tracts_sf %>% pull(tract)) %>%
ca_cvar("tasmax") %>%
ca_period("year") %>%
ca_gcm("ens32max") %>%
ca_scenario("rcp85") %>%
ca_years(start = 2070, end = 2099) %>%
ca_options(spatial_ag = "mean")
kern_trcts_prj_cap
Cal-Adapt API Request
Location(s):
AOI Preset: censustracts
tract(s): 6029000101, 6029000102, 6029000200, ...
Variable(s): tasmax
Temporal aggregration period(s): year
GCM(s): ens32max
Scenario(s): rcp85
Slug(s): NA
Dates: 2070-01-01 to 2099-12-31
Options:
spatial ag: mean
temporal ag (add'l): NA
Fetch values:
kern_trcts_prj_tbl <- kern_trcts_prj_cap %>% ca_getvals_tbl(quiet = TRUE)
## backup: kern_trcts_prj_tbl <- readRDS("data/kern_trcts_prj_tbl.rds")
dim(kern_trcts_prj_tbl)
[1] 4530 8
head(kern_trcts_prj_tbl)
Compute the mean for each tract:
kern_trcts_mean_temp_prj <- kern_trcts_prj_tbl %>%
group_by(tract) %>%
summarise(mean_temp_prj = mean(val))
`summarise()` ungrouping output (override with `.groups` argument)
Now we can join the tables:
dim(kern_tracts_sf)
[1] 151 64
dim(kern_trcts_mean_temp_hist)
[1] 151 2
dim(kern_trcts_mean_temp_prj)
[1] 151 2
kern_tracts_plus_temps_sf <-
kern_tracts_sf %>%
left_join(kern_trcts_mean_temp_hist, by = "tract") %>%
left_join(kern_trcts_mean_temp_prj, by = "tract") %>%
mutate(temp_increase = mean_temp_prj - mean_temp_hist) %>%
select(tract, ciscore, temp_increase)
head(kern_tracts_plus_temps_sf)
Simple feature collection with 6 features and 3 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -13260680 ymin: 4217755 xmax: -13248100 ymax: 4228472
projected CRS: WGS 84 / Pseudo-Mercator
tract ciscore temp_increase geom
1 6029000101 47 5.023285 [K] MULTIPOLYGON (((-13256411 4...
2 6029000102 33 5.034754 [K] MULTIPOLYGON (((-13249411 4...
3 6029000200 50 5.016184 [K] MULTIPOLYGON (((-13253779 4...
4 6029000300 47 5.016184 [K] MULTIPOLYGON (((-13252464 4...
5 6029000400 57 5.016184 [K] MULTIPOLYGON (((-13253750 4...
6 6029000503 15 5.011804 [K] MULTIPOLYGON (((-13260684 4...
Let’s make a choropleth map of the Cumulative Impact score as well as the mean temp_increase:
library(tmap)
tmap_mode("plot")
tmap mode set to plotting
tm_shape(kern_tracts_plus_temps_sf) +
tm_polygons (col = "ciscore",
n = 10,
style = "cont",
palette = "YlOrRd",
colorNA = "grey50",
legend.reverse = TRUE,
title = "CI Score"
) +
tm_layout(main.title = "Kern County Cummulative Impact Score\nCalEnviroScreen",
main.title.size = 0.9,
legend.position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "bottom"))
Mean temp increase:
tm_shape(kern_tracts_plus_temps_sf) +
tm_polygons (col = "temp_increase",
n = 10,
style = "cont",
palette = "YlOrRd",
colorNA = "grey50",
legend.reverse = TRUE,
title = "Mean Temp Increase (K)"
) +
tm_layout(main.title = "Kern County Mean Temp Increases\nHistoric Period - End of Century, 32-ens GCM, RCP85",
main.title.size = 0.9,
legend.position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "bottom"))
NA
NA
To plot CI score and mean temp increase together, we can rescale each one 0..1 and then multiply them together.
kern_tracts_plus_temps_idx_sf <-
kern_tracts_plus_temps_sf %>%
mutate(ciscore_01 = scales::rescale(ciscore),
temp_increase_01 = scales::rescale(as.numeric(temp_increase))) %>%
mutate(csi_temp_idx = ciscore_01 * temp_increase_01)
head(kern_tracts_plus_temps_idx_sf)
Simple feature collection with 6 features and 6 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -13260680 ymin: 4217755 xmax: -13248100 ymax: 4228472
projected CRS: WGS 84 / Pseudo-Mercator
tract ciscore temp_increase geom ciscore_01
1 6029000101 47 5.023285 [K] MULTIPOLYGON (((-13256411 4... 0.6025641
2 6029000102 33 5.034754 [K] MULTIPOLYGON (((-13249411 4... 0.4230769
3 6029000200 50 5.016184 [K] MULTIPOLYGON (((-13253779 4... 0.6410256
4 6029000300 47 5.016184 [K] MULTIPOLYGON (((-13252464 4... 0.6025641
5 6029000400 57 5.016184 [K] MULTIPOLYGON (((-13253750 4... 0.7307692
6 6029000503 15 5.011804 [K] MULTIPOLYGON (((-13260684 4... 0.1923077
temp_increase_01 csi_temp_idx
1 0.1280552 0.07716145
2 0.1439292 0.06089311
3 0.1182271 0.07578663
4 0.1182271 0.07123943
5 0.1182271 0.08639676
6 0.1121651 0.02157021
Plot our index:
tm_shape(kern_tracts_plus_temps_idx_sf) +
tm_polygons (col = "csi_temp_idx",
n = 10,
style = "cont",
palette = "YlOrRd",
colorNA = "grey50",
legend.reverse = TRUE,
title = "CI Score x Temp Increase"
) +
tm_layout(main.title = "Kern County Cummulative Impact Score * End-of-Century Temp Increase \nCalEnviroScreen",
main.title.size = 0.9,
legend.position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "bottom"))
In this example we saw how to join data from Cal-Adapt to other tables using the feature id column. This gives us the ability to combine different types of location-based information for analysis and visualization.